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A real-space renormalization transformation is constructed for lattices of non-identical oscillators 
with dynamics of the general form dcpk/dt = uik + fik{4>h4'k)- The transformation acts on 
ensembles of such lattices. Critical properties corresponding to a second order phase transition 
towards macroscopic synchronization are deduced. The analysis is potentially exact, but relies 
in part on unproven assumptions. Numerically, second order phase transitions with the predicted 
properties are observed as g increases in two structurally different, two-dimensional oscillator models. 
One model has smooth coupling fik{4'i,4'k) = '4>{4>i~4>k), where ifi{x) is non-odd. The other model is 
pulse-coupled, with fik{4>i, 4>k) = S{4>i)(pi4>k)- Lower bounds for the critical dimensions for different 
types of coupling are obtained. For non-odd coupling, macroscopic synchronization cannot be ruled 
out for any dimension D > 1, whereas in the case of odd coupling, the well-known result that it can 
be ruled out for D < 3 is regained. 

PACS numbers: 05.70.Fh, 05.10.Cc, 05.45.Xt 
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I. INTRODUCTION 

The study of synchronization in large oscillator net- 
works has been a thriving field of research ever since the 
classic work by Winfrce in 1967 [l,]. Even so, there are 
still basic questions that await satisfactory answers. One 
such question is when and how macroscopic synchroniza- 
tion occurs in lattices of non-identical oscillators. This is 
the subject of the present paper. 

Part of the charm of the study of synchronization is 
that the applications are very diverse [2;, 3, 4] . Rhythmic 
activities in living organisms are, in many cases, gen- 
erated by the collective oscillation of a large, synchro- 
nized assembly of pacemaker cells. Examples include the 
beating of the heart [s], locomotion the circadian 
rhythm 0] , and the peristaltis of the small intestine Q . 
There is also growing evidence that large-scale synchro- 
nization among neurons is crucial in the interpretation of 
sensory data and in conscious perception 9]. Epileptic 
seizures correspond to an abnormal degree of synchro- 
nization [lo{. On a larger scale, synchronization can be 
seen in groups of organisms. Swarms of fireflies may flash 
in unison the chirping of crickets in a field waxes 
and wanes in partial synchrony [T^ . an audience may 
spontaneously start to clap in unison [l^, females liv- 
ing together synchronize their menstrual cycles p^ . It 
has also been realized that synchronization is an essen- 
tial concept in the dynamics of spatially extended animal 
populations fis']. Examples from outside biology include 
synchronization in power grids among lasers 16], os- 
cillatory chemical reactions and in arrays of Joseph- 
son junctions p^. 

In most applications, there will inevitably be some 
variation among the oscillators, for instance in the natu- 
ral frequency with which they oscillate when isolated. In 
this situation, macroscopic synchronization means that 
the order parameter r becomes non-zero, where 



Here, M is the size of the largest group of oscillators that 
attain the same mean frequency and N is the total num- 
ber of oscillators. If the network has spatial structure, 
the M synchronized oscillators typically form a percolat- 
ing cluster [l^, |23]. The mean frequency Qk of oscillator 
k is defined as 



flk = lim 4>k{t)/t, 

t — >oo 



(2) 



where (j>k is the phase of k. The existence of the above 
limits has to be assumed ''2^. 

In theoretical work, the description of each oscillator 
must be simple to enable the study of large networks. 
Kuramoto and others 22] introduced the so called phase 
reduction technique, and showed that in the limits of 
small coupling between oscillators and small variation of 
natural frequencies, the phase (jjk is sufficient to describe 
the state of each oscillator fc, and the network dynamics 
is given by 



dt 



N 

1=1 



(3) 



The constant is the natural frequency, g is the cou- 
pling strength, and ipik{x) is a 1-periodic function. An- 
other situation where a phase description is sufficient is 
in the limits of short, pulse-like interactions and strong 
dissipation. The quick reduction of phase space volumes 
then ensures that after one perturbation from a nearby 
oscillator /, oscillator k returns close to its limit cycle 
before the next perturbation, and we may write 



d(j)k 
dt 



N 



—— ^Ljfc-t-g > (5[mod(0;,l)]v?ifc((?!)fe), 



1=1 



(4) 



r = lim M/N. 



N 



(1) 



where 5{x) is the Dirac delta function. This is often a 
good description in biological applications, where the in- 
teractions, for instance, may consist of electric discharges 
or light flashes. If we define 0^ as a cyclic variable, 
(/>fc S [0, 1), we can replace 5[a\od{(f)i , 1)] with 5{(t)i). The 
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1-periodic function ipik{x) is often called the phase re- 
sponse curve. 

Transitions to macroscopic synchronization are simi- 
lar to phase transitions in equilibrium systems. The 
two main methods to analyze phase transitions are to 
make a mean-field description or to use a renormaliza- 
tion group. So far, most attempts to gain understanding 
of macroscopic synchronization among non-identical os- 
cillators have assumed that the oscillators are coupled 
all-to-all. This is the mean-field description. To use a 
renormalization group, on the other hand, is the natural 
way to gain understanding of phase transitions in lattices. 

This is the method used in this study. A real-space 
renormalization scheme is developed that is potentially 
exact. It is tested numerically on two structurally differ- 
ent models. At its present stage of development, however, 
the scheme has to be called heuristic since it relies in part 
on unproven assumptions. 

Before describing the approach, let me review very 
briefly the current state of knowledge about transitions 
to macroscopic synchronization in mean-field and lattice 
models. 



II. REVIEW OF RELATED WORK 



Mean-field models 



As a special case of Eq. ([3]) , Kuramoto introduced 
the mean-field model 



mean-field model, r and R typically becomes non-zero at 
the same critical coupling gc- In a lattice model, waves 
in the phase field may be expected [E [3, i, II [H, [M HI 
even if the frequencies are synchronized, so that r > 
but i? = 0. 

Ariaratnam and Strogatz [27] studied Winfree's origi- 
nal model 



dt 



1=1 



N 



(8) 



in the special case d{x) = 1 + cos(27ra;) and (p{x) — 
— sin(27ra;) . This model is similar to the model ^ , with 
the smooth 1-periodic influence function 9{x) replacing 
the delta pulse. The authors were able to obtain the 
phase diagram in the plane spanned by g and 7, where 
is uniform with support [1 — 7, 1 -I-7]. Apart from the 
phases with r — R — Q and r = R = 1 , there is a phase of 
partial synchrony with r < 1 and i? < 1, and also phases 
with partial or complete oscillator death {d<j)kldt = 0). 

Tsubo et al. [28] studied a similar model, but let 
the disorder reside in the phase response curves ipik{x), 
whereas the natural frequencies were identical. With 
Vik{x) = cos(7rafe) — cos(27ra; — irau), where ak is a ran- 
dom number from a uniform distribution with support 
[oniin, Q-max], they found a discontinuous transition to 
macroscopic synchronization in the phase plane spanned 
by ctniin and amax7 in contrast to the continuous transition 
in the Kuramoto model, as expressed in Eq. ([7]). 



N 



dt 



= + ]V ^ sin[27r((/)/ - (/)fc)]. 



(5) 



1=1 



because of its analytical tractability. Instead of r, Ku- 
ramoto studied the order parameter 



R = lim lim 



N 
k=l 



(6) 



and found that there is a critical coupling strength gc 
such that 



i? = 

R(x{g 



gc 



9 < 9c 
1'/' 9>9c 



(7) 



If each LOk is chosen independently from a density func- 
tion that is unimodal and symmetric about its mean 
fi, the critical coupling is given by gc = 2/7rI?^(/i). 

Since the original work b y K uramoto, the analysis of 
model ([5]) has been refined [2^. Also, it turns out that 
the exponent 1/2 in Eq. ([7]) changes to 1 as soon as 
non-odd harmonics are added to the coupling function 
(fik{x) = sin(27rx) [2^. 

Note that the order parameter R measures the degree 
of phase synchronization, whereas r [Eq. ([1])] measures 
the degree of frequency synchronization. A non-zero R 
implies a non-zero r, but the opposite is not true. In a 



B. Lattice models 

The analysis of oscillator lattices is harder than that of 
mean-field models, and less progress has been made. For 
cubic lattices with dimension D and dynamics of form ([3]) 
with f>ik{x) — ip{x) and odd coupling, ip{—x) = —ip{x), 
Daido ruled out states with r > for £> < 2 [291]. 
Daido obtained this result using renormalization-like ar- 
guments. With similar methods, Strogatz and MiroUo 
[30| were able to prove that whenever D^^ has non-zero 
variance, states with r = 1 are ruled out for any finite 
D. In addition, states with < r < 1 cannot have syn- 
chronized clusters which contain macroscopic cubes (with 
volume V — aN, < a < 1). Thus, for odd coupling, 
macroscopic synchronization may occur only if D > 3 
and can only be partial, with sponge-like synchronized 
clusters. Whether such states actually exist is still an 
open question. The numerical evidence is inconclusive in 
my view [H [H, [HI, [13 . 

Kopell and Ermentrout were the first to point out that 
non-odd coupling facilitates synchronization fssj . For an 
oscillator chain [D — 1) of form ^ with <fik{x) — fix), I 
studied the case when has finite support [o^min, t^^max]- 
For models with ip{Q) = and (p'{0) > 0, there is then 
a critical coupling gc at which a discontinuous transition 
from r = to r = 1 takes place [sj. I found that gc = 
('j-'max — '^min)/ \d{x)\, whcrc the denominator is a 
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mesure of the "non-oddity" of </5(a;), vanishing for odd 
couphng such as (p{x) = sin(27ra:). A similar result was 
provided for a model of the form ([4]), which can be seen 
as inherently non-odd due to the sequential interaction 
of two oscillators via pulses [s^ . 

Since macroscopic synchronization is possible for D = 
1 for non-odd coupling, it is expected to be possible for 
all D > 1. However, no proofs have been obtained, to 
my best knowledge. For D = 2, me and my co-workers 
[20| offered numerical evidence for a continuous, second 
order phase transition to r > in a model of form 

In equilibrium systems, there is typically an upper crit- 
ical dimension, above which a lattice model shows mean 
field critical behavior. Hong, Park and co-workers (3^ 
have re-examined the lattice version of the Kuramoto 
model ([5]), and claim that _D = 4 is the upper critical di- 
mension, above which critical exponents take mean field 
values and macroscopic frequency- and phase synchro- 
nization appear at the same critical coupling. However, 
the results by Strogatz and MiroUo [30[ indicate that the 
upper critical dimension is infinity for this model, since 
they ruled out states with r = 1 for any finite D and 
any with non-zero variance, whereas such states exist 
in the mean-field model ^ when has non-zero vari- 
ance, but finite support. It is conceivable that the phase 
transition structure of oscillator networks is richer than 
in equilibrium systems, and cannot be fully captured by 
the concepts used there. 

Lattice models of oscillator networks are closely re- 
lated to spatial continuum models. This is the natural 
way to describe the oscillatory Belousov-Zhabotinsky re- 
action [sH and the smooth muscle tissue in the intestine 
Q. It may also be an adequate model of a large piece 
of oscillatory cardiac muscle, even though it consists of 
discrete cells. The preferable mathematical description 
is given by the Ginzburg-Landau equation (GLE) (3l.[2^. 
where the state at each point in space is given by a com- 
plex number, encoding both the phase and amplitude of 
oscillation. The GLE corresponds to a lattice of iden- 
tical oscillators. Using a field theoretic renormalization 
group, Risler and co-workers have performed a thorough 
analysis of synhronization transitions in the GLE with 
noise [37| . The noise is assumed to be uncorrelated in 
space and time. In contrast, random natural frequencies 
correspond to "noise" that is uncorrelated in space, but 
quenched in time. This makes the problem much more 
difficult in the continuum formulation. In particular, dis- 
continuities arise in the phase field whenever frequency 
synchronization is not perfect (r < 1). 



III. MODELS AND METHODS 

In the analysis, models of the following form are con- 
sidered: 



Model 1 



Model 2 



i.5r 
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dt 



^fc+.9l] /ifc(0/,0fc), fc = l,...,iV. (9) 



1 



FIG. 1: Coupling functions in the two test models. Model 1 
is of form ((3} with ipik{x) given by Eq. (|12|l . Model 2 has 
form Q with ipik{x) given by Eq. ([13} 



Here, (j)k € M is the phase of oscillator fc, ujk is its natural 
frequency, g is the coupling strength, and rik is the set of 
fc's nearest neighbors. The analysis is restricted to cubic 
lattices of dimension D. The coupling functions fik are 
assumed to be 1-periodic in each argument. With this 
restriction, the phases (j)k are allowed to grow linearly to 
be able to count the number of cycles, that is, the largest 
integer smaller than (j)k{t) — 0fc(O). Since no further as- 
sumptions are made, the results are expected to apply 
(at least) to all models of this form. All the coupling 
functions in the models referred to above have the form 
given in Eq. ©. 

Let us define the ensemble 

£ ^£{g,V^,Vf,V^^o),D,N) (10) 

of realizations of systems ([9]), where ujk are indepen- 
dent random numbers from the density function T?^, each 
fik is chosen from Vf, and the initial condition 0(0) = 
[(/)i(0), . . . , (^Ar(O)] is chosen from I?0(o). Quenched disor- 
der is introduced by and Dp 

To give the coupling strength g a clear meaning, Vf 
should be chosen so that 

(/ / |//fc(0i,0fe)M0id'^fe)D, =1, (11) 

Jo Jo 

or so that it fulfils a similar condition. Alternatively, one 
may drop g as an argument of £. 

To test the theoretical predictions, numerical simula- 
tions of two specific models with D = 2 are performed. 
The first model has the form ^ with 

Vik{x) — sin(27rx) -f- ^ sin^(27rx), (12) 

(Fig. [1]). This model will be referred to as Model 1. The 
density function T^i^ is uniform with support [1.0, 1.5] and 
I?0(o) is uniform in the interval [0, 1]. Forward Euler inte- 
gration is used, with At = 0.05. The motivation for this 
model is that it is similar to the Kuramoto model, but 
has a non-odd term to allow macroscopic synchronization 
for D = 2 (see below). 

The second model has the form ^ with 

-X, X < 0.4 
'fik[x) = '{ 9a; -4, 0.4 < a; < 0.5 , (13) 
1-a;, 0.5 < a; < 1 
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FIG. 2; A block-oscillator transformation (|14|l with scale fac- 
tor 6 = 2 and implicit change of length scale r' — r/b. 



(Fig. [T]). This model will be referred to as Model 2. V^-i 
is uniform with support [1.0, 1.5] and ^'^(o) is uniform in 
the interval [0, 1]. The same integration method as in 
Refs. [23| and Q is used. The piecewise linear phase 
response curve expressed by Eq. (fT3|) has the same bipo- 
lar characteristics as the curve 1^3(2;) = — sin(27ra;) used 
by Ariaratnam and Strogatz [27]. This type of response 
to external perturbation is found in many biological ap- 
plications [3]. 

The details are arbitrary, but these models are chosen 
since they have been studied previously, model (|12p for 
D = lm Ref. [34] and model m for D = 1 in Ref. [H 
and in the case = 2 in Refs. [2C|] and (24j . 

Unless otherwise stated, simulations of Model 1 are 
carried out with lattice size 300 x 300, whereas for Model 
2 the lattice size 500 x 500 is used. The larger lattice 
size is needed to resolve phase 2 in Model 2 (see below). 
Periodic boundary conditions are used, unless otherwise 
stated. In and around critical regions, a transient time 
oi t = 100000 — 200000 is used before measurements are 
done. For g well below critical values, shorter transient 
times are used. Mean frequencies flk are approximated 
by taking the mean of d^fc /dt during a time interval At = 
1000 after the initial transient [2l[. 

To identify frequency clusters, the lattice is scanned. 
An oscillator k is considered to belong to the same 
cluster as a previously scanned neighbor oscillator I if 
]r2; — fifc] < 0.001. If two such neighbor oscillators I and 
I' are preliminarily judged to belong to different clusters 
C and C", but both fulfil the above inequality, C and C" 
are identified as two parts of the same cluster. 



IV. THEORETICAL APPROACH 

The first step in the renormalization scheme is to define 
a block-oscillator transformation pi, : M'' +^ ^ (Fig. 



(14) 



We get a coarse-grained version of the lattice and inter- 
pret the phase (f)j as the state of block oscillator j. Ap- 
plying pb to all j, we may define the scale transformation 



A 

V 
S 



{W)} 

A 
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FIG. 3: Relation between the renormalization transformation 
Rt and the block-oscillator transformation pt- An ensemble 
£ produces an infinite set {(pit)} of time series (j}{t). Only 
if such a set {^(t)} inversely determines £ (dashed vertical 
arrows), is £' uniquely determined by {^(t')}. Thus Rt does 
only exist in this case. 



Ph 



as 



m'),t']^Ptm),t], 



(15) 



with (j) = (01, . . .,(j)N) and (i> = (0i, . . . ,0jv/6")- 

Let us discuss these transformations in a bit more de- 
tail. To be able to interpret Pf, as a scale transformation, 
note that it must fulfil the group property 

Pb.Pb, = Pb,b,- (16) 

Regarding pb, I restrict the interst to linear transforma- 
tions of the form 



M{b, D) 



( 4>k, \ 
V t J 



(17) 



where 



and M{b, £)) is a 2 X (6^ 
M{b, D) = 



•k^o is some list of the phases in block j, 
l)-matrix with 



A B 
C 



(Here, A, B, and C are functions of b and D.) Lin- 
ear transformations with Mi_i = . . . = Mi b^ — A are 
considered since 4>j is intended to be a kind of arith- 
metic mean of the phases 0fc in block where all 4>k are 
treated in the same way. Such a choice of (pj is justi- 
fied if the phases are interpreted to be linear variables, 
instead of cyclic ones. There is no reason to let the trans- 
formed time t' depend on the phases, and therefore I set 
Maa = . . . = Ah.bo = 0. 

Next, let {H[(f){t)\) ^ be the ensemble mean of i?[0(i)], 
where H \s a. functional of the time series (/>(i). We then 
seek an ensemble 

E' = £[g', Vf, , V^, , I?^,(o) , D, N/b"") (18) 

[c.f. Eq. (dni)] such that 



(//[</.' (t')])£' - {H[m])e 



(19) 



for any functional H. In words, all statistical quantities 
produced by the desired ensemble £' should be the same 
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FIG. 4: Projection of {5} to the plane spanned by coupling 
strengh g and the variance cr^ of natural frequencies. The 
transformation [Eq. (|14p ] is chosen so that a critical fixed 
point E* with finite g* and (o"5)* may appear. The flow under 
Rh is intended to be such that ensembles on the critical line 
Sc are attracted to f *, and that cr^ is invariant at jr = 0. 
Numerically, two ensemble families Eg are studied, Model 1 
and Model 2. Each of these seemingly becomes critical at 
some coupling g = gd- 



as those given by the transformed phases {^(t')]^ which 
in turn are determined by E. The relation 



(20) 



defines the renormalization transformation i?;,, assuming 
that E' exists (Fig. EJ. 

Naively, instead of working at the ensemble level, we 
could have looked for an evolution equation for the trans- 
formed phases ^. If we would have been successful, we 
could have written d^jjdt' ~ ^j+g'J2ien.i flji'Pi^'I'j) 
in compressed form, d(j)j/dt' = ujj+h'{(j))- However, since 
Pb is not invertible, we have to express the interaction in 
the original phases, i.e. h' = h' {(/)), and we do not get an 
evolution equation of the transformed phases in closed 
form such as Eq. 

Let us write 



bk/dt 
b'Jdt' 



hki4>) 

Km 



(21) 



for the original and transformed ensembles E and f , re- 
spectively. In the following, I also adopt the notation 
E[a;] = {x)s and let Var[a;] and Cov[a;, y] be the ensemble 
variance and covariance, respectively. 

To be able to extract information about critical behav- 
ior from i?f,, the transformation pi, has to be chosen so 
that there may appear a non-trivial fixed point ensemble 



E* = RbE* 



(22) 



in the limit N ^ oo (Fig. [5]). In this study I look for, 
and assume the existence of, a fixed point E* for which 



the variance of natural frequencies and the variance of 
the interaction exist, that is 



< oo 



(23) 



< Var*[/ifc] < oo ' 

A finite Var* \hk\ implies a finite fixed point coupling 
strength g* |38| . 

Further, E* should attract an ensemble Eg^-^ belonging 
to a family Eg that passes a transition to macroscopic 
synchronization at the critical coupling g = gd- The 
behavior of Eg^ will then be the same as that of £* at 
large scales and after long times. 

With this in mind, pb is chosen to fulfil three condi- 
tions, in addition to Eqs. and P?|) . Before stating 
these conditions, let me introduce a few quantities. 

First, let m{t) and moo be mean attained frequencies: 



m{t) 



{dcpk/dt)£ 
limf^oo ■m{t). 



(24) 



The limit limf^oo rn{t) exists at the presumed fixed point 
E* according to assumption ((23| . In fact, it follows from 
Eqs. (PTjl and that the two first moments of the 
distribution of attained mean frequencies fifc exist at E* 
[Illli: 



E*[Q.k\ = (moo)* < oo 
< Var*[r2fc] < oo. 



(25) 



Further, let k be the mean wave number. Since the 
phases are allowed to be linear variables, the wave nature 
of the phase landscape in a given lattice £ may only be 
manifest if a suitable integer is added or subtracted to 
each Writing Q = (gi, . . . , g^r) and _ ^ _j_ j 
define 



K = lim Min 

t — ''OC 



(26) 



In other words, Q should be chosen so that the mean 
phase difference between neighbor oscillators is mini- 
mized, and this phase difference is k. The lattice mean 
(. . .) c is expected to become equivalent to an ensemble 
mean (. . .)£ in the limit TV — > oo. (Otherwise an ensem- 
ble mean can be added in the definition.) 

Let m and h be the corresponding mean frequency and 
wave number in the transformed lattice. The three con- 
ditions that guide the choice of pf,, apart from Eqs. (fTBl) 
and (|17p. are then: 

1. There is a (non-empty) set of ensembles Si, such 
that if £" € El, then moo — moo for any h. 

2. There is a set of ensembles E2 C Ei, such that if 
£ e S2, then k — k for any b. 

3. There is a set of ensembles E with no coupling (g — 
0) such that if cr^ is finite and non-zero, then <t^, ~ 
cr? is finite and non-zero for any b. 
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il) < b = 2 ^ 
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I J 

FIG. 5: Illustration for the argument why condition 3 is ful- 
filled for the block-oscillator transformation pb [Eq. (|28|l ]. A 
wave is moving in the positive a;-direction in the case D = 1. 
Two neighbor block-oscillators i and j with size 6 = 2 are indi- 
cated. Corresponding individual oscillators in the two blocks 
are paired as (^i. , ) and (02^ , 02^ ). The oscillators in each 
pair are placed at distance b = 2 from each other. 



It is clear that the two first conditions are necessary. The 
critical fixed point £* I hope to construct belongs to S2. 

To see why condition 3 is necessary, look at Fig. 
m Assume first that the condition is broken and that 
limb_+oo cr^, — for all ensembles with no coupling. Then 
the unstable manifold U oi £* bends down to the origin. 
If Sc would still be a stable manifold of £*, closed flow 
lines would appear, which is impossible since correlation 
lengths are always reduced a factor b each time Rb is ap- 
plied. Thus the flow along the critical line Sc changes 
direction, and the critical properties at g^i are no longer 
given by those of £* . In other words, £* becomes irrele- 
vant. 

Assume instead that limt_>oo f^' = 00 at g = for 
all ensembles with no coupling. I will give a plausibility 
argument why this is not consistent with the existence of 
a fixed point £* of the desired kind. From Eqs. (|17p and 
m, d4>,/dt' = iA/C)j:,^^Lu, + {A/C)j:,^^hk + B/C. 
Since g = corresponds to /i/j = 0, we have d(t)j/dt' — 
i^j)g=o + (^/C)Efeej^fc- Here, (wj)g^o is the trans- 
formed natural frequency of block j in the ensemble ob- 
tained when g is replaced by in the original ensem- 
ble £{g, ■ ■ .). Taking the ensemble variance and apply- 
ing the equation to the presumed fixed point £*{g*, . . .), 
we may write Va-r* [d(pj /dt'] — (cr^/)g_Q + R, where I 
do not specify the rest term R for clarity. The left 
hand side of this equation must be finite for any 6, 
since Yar[d(j)j / dt'] = VaT[d(j)k / dt] at a fixed point and 
Yai* [d(f>k/dt] exists [S^. This condition would be hard 
to fulfil if (cr^,,) „ ^ 00 as & — *■ 00. Then the term R 
must sensitively balance this divergence. 

Straightforward algebra shows that the only block- 
oscillator transformation of form (fTTl) that has the 
group property and satisfies conditions 1 - 3 is the 
one with 

A = 6-^-1 

B = - 6-i)moo (27) 

c = 6-^/2-1, 



or, explicitly, 

r ~4>,{t') = 6-^-1 Efce, Mt) - b-\l - h-^l-')m^t 
\t' = b-D/^-H 

(28) 

(In fact, I only demonstrate that this transformation sat- 
isfies condition 2 in a restricted sense, to be described 
below.) 

Note that if we are interested in critical ensembles for 
which cr5 does not exist, the three conditions, and hence 
transformation . ^hould be modified. This point is 
discussed by Daido [29i|. This in turn affects the critical 
properties and the critical dimension. I do not deal with 
these cases explicitly in this paper. 

Instead of proving the uniqueness of transformation 
, let us examine to what extent it fulfils conditions 1 

- 3. 

Regarding condition 1, we get 

rhoo = TOoo (29) 

for any b and any £ by direct evaluation of 
\imt^oc.{d^j/dt')£, using Eq. 

Condition 2 is fulfilled in the following restricted sense. 
If we may choose Q [Eq. ((26| ] such that a plane wave 
moving along a principal axis appears in the phase field 
(jjiQ) (^j-'j ^ then its wave number will reamin the same in 
the transformed field 0('^^(r'). Assume that i and j are 
two neighbor block-oscillators along the relevant princi- 
pal axis (Fig. [5]). Then, dropping the superscript (Q) for 
brevity, 

b° 

0,-^, = 6-^-1 ^(0fc^.-0fcj, (30) 
fc=i 

where ki and kj are the oscillators at corresponding posi- 
tions in block i and j, respectively. The distance between 
these is b, and thus (c^fc^. — (l)ki)£. = b{4>i — 4'k)c, where 
k and / are neighbor oscillators along the principal axis. 
Consequently, taking the lattice mean of Eq. pO)) . we 

get - cpi)c = (0i - <^fc)£- 

Turning to condition 3, let us use Eqs. (p8)l and ([29]) 
to write 

f---'-°'^i:<t--)^ ,31, 

kej 

At 5 = we have d(j)j/dt' — ujj and d(j)k/dt = ojk- Also, 
Woo = rnoa — fl = ^1 hy Eq. so that ujj — jl — 

6-^/2 Y,^^^[ujk - fi). Thus al = al for all b. Condition 
3 is then fulfilled since Vi^i — at g = 0. In the limit 
b —^ 00, V^i becomes a Gaussian by the central limit 
theorem, with the same mean and variance as Vi^. 

Transformation (pS)) is the only block-oscillator trans- 
formation of form p7|) that enables a non-trivial fixed 
point £* of the desired kind. There may be acceptable 
transformations that do not have form ([T7]). This is not 
essential. Fixed point properties derived from any pt that 
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gives rise to a fixed point £* of the desired kind have to re- 
flect critical properties of model © , if the corresponding 
family of ensembles £g pass through the critical surface 
Sc as couphng strength g is varied (Fig. U]). 

To gain some information about Rb^ let us try to ex- 
press 



d(f)j/dt' — ojj + hj{4>), 



(32) 



Model 1 



0.18 




9c2 Model 2 
0.161 



N 



0.30 



9ci 



0.60 



where ujj is a constant that can be interpreted as the 
natural frequency of block-oscillator j in the transformed 
lattice, and hj{(f) can be seen as the interaction term. To 
allow such an interpretation, hj{(j)) has to be zero when 
block j is decoupled from the rest of the lattice. Let 

m., = lim b-°y^[d(l)k/dt]^, (33) 



FIG. 6: The portion M-/N of the lattice occupied by the next 
largest frequency cluster as a function of coupling strength 
g. The maximum of M-/N corresponds to gci and it drops 
almost to zero at gc2- For Model 1 it is found that Qci ~ 0.23 
and gc2 ~ 0.28, and for Model 2 g^i ~ 0.50 and (?c2 ~ 0.56. 
Average and standard deviation of 7 realizations of Model 1 
and 3 realizations of Model 2 are shown for each g. 



where is the value of x when block-oscillator j is de- 
coupled from the surroundings. Further, let an under- 
lined variable denote an intitial condition mean: 

x = {x)v^^oy (34) 

In a sense, m^- is then the natural frequency of block- 
oscillator j in the original lattice. Using Eqs. ([9]) and 
(PS)) . we may express 

b-^'/'Eke, [g [Eien, - (m, -^fe)] • 

(35) 

Let us interpret: 

Uj = TOoo + b^^^irUj - moo), (36) 

and 

/i,(<^) = 6-^/2 Efce, {9 [Eien, fiki<l>i,M] - (m, - ^fe)} 

(37) 

Here, /^"^ is the initial condition mean of the coupling 
function fik{4'ii4'k) as t — > oo and block j is decoupled, 
and — g^ifzn^ /;"^.- The interaction hj in Eq. (157)) is 
not strictly zero when j is decoupled from other blocks. 
However, the initial condition mean is zero, in the 
limit t ^ oo. Thus, these identifications can be used to 
deduce asymptotic critical behavior of initial condition 
averaged variables, but nothing else. 
I formulate the following conjecture: 

{HWM')]}e' = {Hih^ims ^^^^ 

for any functional _ff as t — > oo. A number of critical 
properties follow from Eqs. (HH), (gS]), ((351), and the fixed 
point condition £' = £ = £*. 



V. RESULTS 
A. Phase diagrams 

I want to compare each theoretical prediction with nu- 
merical results from the two test models. Model 1 [Eq. 
(HH)] and Model 2 [Eq. in the case D = 2. To do so, 

it has to be demonstrated that there are critical points 
in these models, and the critical coupling strengths gd 
have to be identified. 

Previously, two critical couplings gd and gc2 were 
found in Model 2 [23|- The system seemingly becomes 
critical at g ~ gd and almost perfect synchronization 
settles a.t g — gc2- There are still isolated oscillators that 
never fire for g > gc2, and thus they are not synchronized 
to the rest of the lattice [l^l- The two critical couplings 
separate three phases, phase 1 {0 < g < gd), phase 2 
(ffci < 9 < .9c2), and phase 3 {g > gc2)- These con- 
clusions were reached by looking at the distribution of 
cluster sizes. At g = gd, this distribution seems to obey 
a power law. This is true in phase 2 also, if the macro- 
scopic cluster is disregarded. By estimating gd and gc2 
for different values of A^, it was argued that the three 
phases are not a finite size effect, but persist as TV — > cx). 

Simulations with Model 1 suggest that it behaves qual- 
itatively in the same way. Instead of showing cluster 
size distributions, we look in Fig. [S] at the quantity 
M_/7V, where Af_ is the size of the next largest clus- 
ter. M^/N is expected to peak at gd- Above gd the 
largest, percolating cluster grows in size as g increases 
further, whereas the other clusters become smaller. At 
5c2, M-/N should drop close to zero. In this way, it is 
estimated that gd ~ 0.23 and gc2 ~ 0.28 for Model 1, 
whereas gd ~ 0.50 and gc2 ~ 0.56 for Model 2 (Fig. 0. 
Figure [7] shows frequency landscapes ri(r) in each of the 
three phases (2l| . 

For both models, phase 2 is rather narrow, but clearly 
distinguishable. To establish the existence of phase 2 
even more clearly, complementary simulations were made 
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(0 
Q. 



Q. 




9 = 0.20 (1.265,1.310) g = 0.45 (0.870,0.910) 





g = 0.23 (1.275,1.300) g = 0.50 (0.900,0.925) 




0) 



CO 
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g = 0.30 (1.280,1.300) g = 0.60 (0.932,0.942) 



FIG. 7: Frequency landscapes Q{r). Frequency is coded ac- 
cording to the bracket (tJmin i^max)- An oscillator k is colored 
black if Qk is less than Wmin and white if it is higher than 



o.sr 



Model 1 g,;-] 



9c2 



Model 2 



0.8 



Phase 1 



















Phase 1 










Phase 3 



0.4 



0.8 



FIG. 8: Phase diagrams. The same quantity M~ /N as in Fig. 
[6] is used to identify Qd and gc2- For model 1, T>lj is uniform 
with support [1, 1 + 7]. For model 2, V^-i is uniform with 
the same support. Phase 1 is subcritical with microscopic 
frequency clusters only. In phase 2, there is one macroscopic 
cluster. In phase 3, almost all oscillators synchronize their 
freqencies. All three phases seem to persist even if T>uj has 
tails (see text). For Model 2, it is impossible to resolve phase 
2 in the data obtained with 7 = 0.3. I cannot decide whether 
phase 2 extends down to the origin, or if there is a triple point. 



for wider and narrower (Fig. [8]). 

The hypothesis that will be tested is that the curve 
(jfci(o'S) is identical to the critical curve Sc, which is also 
the stable manifold of a critical fixed point 8* (Fig. [5]). 
Thus I identify — gd ■ It is a delicatequestion whether 
the entire phase 2 is critical. In Ref. [20] I hypothesized 
that this is so, based on the cluster size distribution and 
the temporal instability of clusters, even after very long 
times. The large sample-to-sample fluctuations seen in 
Fig. [5] further strenghten this idea. The matter is dis- 
cussed further below. 

To test whether phase 3 exists even if D^i has tails. 
Model 1 is simulated with Gaussian natural frequencies, 
with mean fj, — and variance a^, — 1/48, i.e. 



2?^ =AA(0,l/48). 



(39) 



The variance is chosen to be equal to the variance of the 
original, uniform Vi^. I estimate gd — 0.23 and gc2 — 
0.28. Phase 3 is entered even if the oscillators with the 
most extreme natural frequencies do not synchronize to 
the rest of the lattice. Model 2 is simulated with the 
Rayleigh density function 



V,. 



47r(w~ 
0, 



l)e-27r(a;-i-l) = 



> 1 
-1 < 1. 



(40) 



In this model, it is found that gd ~ 0.55 and gc2 ~ 0.75. 
Phase 3 is entered even if the slowest oscillators do not 
synchronize to the rest of the lattice. 

I hypothesize that the transition to phase 3 is discon- 
tinuous, since the distribution of cluster sizes seems to 
collapse discontiuously at gc2- However, more detailed 
studies are needed to establish the nature of this transi- 
tion, and to be able to define phase 3 precisely. 
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FIG. 9: The frequency correlation length [Eq. (|44|l ]. Ac- 
cording to Eq. (|43p . is expected to drop to zero at g — g^i. 
The same number of realizations as in Fig. [6] are used. Since 
linear interpolation is used, zero correlation between neigh- 
bors gives — 1 — (dotted horizontal line). 



B. Frequency correlations 



We have £[17'^] 
and (129]) and 



E[fifc] from Eqs. dUD 



Var[f];] = Var[Ofc] -h 6--° J2 Cov[nk,i^k'] (41) 

k,k'£j, k'=ik 

from Eq. (I3ip '^T]. At a fixed point £*, the sum has to be 
zero for any b if Var*[rifc] < oo, which is the case treated 
here [Eq. (^5)) ]. Therefore we must have Cov*[Qk, ^k'] = 
for all fc 7^ fc. Introducing the pair correlation function 

ro(r) = Cov[nk,nk']ik'-k\=r/yaT[nk], (42) 

it is concluded that 



TUr) = 0,r>l. 



(43) 



Note that even if r^(r) = 0, the fifcS do not have to be 
independent at a critical fixed point. Rather, clusters of 
oscillators which run at the same frequency are expected 

(Fig. ED. 

Figure [5] shows the correlation length ^n, defined by 
the relation 



rn(fn) 



(44) 



(Note that this is not the standard way to define a corre- 
lation length, thus the hat symbol. Normally it is defined 
as the rate of exponential fall-off at large r. C.f. Eq ([82]) . 
The quantity is used here since it turned out to be 
more stable, given the fluctuations of rn(r).) In Model 

1, drops significantly just above the estimated value of 
gd, with comparably small sample variations. In phase 

2, seemingly increases again, even if large variations 
make such a conclusion uncertain. In Model 2, falls 
steeply towards zero just below gd, and stay very close 
to zero for g > g^i with very small variations. 

Thus, Model 2 supports the theory better than Model 
1 does. However, the drop of in Model 1 may occur 
exactly at gd, given the uncertainty in its estimation due 
to the large sample variations of the cluster sizes (Fig. O. 

No indications of negative correlations have been ob- 
served. 
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FIG. 10: Density functions Dn of attained frequencies Slfe. in 
Model 1 and Model 2. The vertical dashed lines in panels 
(a) show the mean frequency E[r2fc]. Panels (b) show the low 
frequency tails. The dashed lines correspond to a relation 
'Dn{Qk) oc (fipoak — fifc)"*. The critical coupling strengths 
are g^i = 0.23 for Model 1 and gd = 0.50 for Model 2. 



C. Frequency distribution 

Consider the density function of attained frequen- 
cies flk [131 • I do not attempt to deduce (Vq)* but dis- 
cuss some of its basic properties. 

We have already concluded that assumption P5)) im- 
plies that E*[rife] and Var*[rife] both exist. Put differ- 
ently, if the critical properties of Model 1 and Model 2 
are to be the same as those deduced for the fixed point 
£*, then E[rii:] and Var[r2fe] should exist at g = gd- In- 
finite Var[r2fc] at g = gd and finite Var*[rife] make neces- 
sary negative frequency correlations [Eqs. (|41]) and (j42|) ] 
along the critical line Sc (FigH]), which have not been 
observed in simulations. 

At g — gd, the order parameter r [Eq. (jT])] becomes 
non-zero, that is, a finite portion of the oscillators at- 
tain identical frequencies, so that an infinitely high spike 
develops in I?^ as g approaches gd from below: 



Max[I?o] 
Let us use Eq. (1^ to write 



lim 

-(9cl). 



OO. 



(45) 



= E[r!fe] + 6-^/2 Y,i^k - mk])- (46) 
kej 

At a critical fixed point, the emerging spike should not 
move when pt is applied, so that 



(f^pcak)* — E*[ili 



(47) 
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FIG. 11: Density functions V^, of attained renormalized fre- 
quencies Clj. Equation (|46|l is used for the numerical renor- 
malization with different scale factors b. In phase 1 {g — 0.20 
for Model 1 and g = 0.45 for Models 2), D^i approaches a nor- 
mal distribution correponding to a trivial fixed point at g = 0. 
This does not seem to be the case in phase 2 {g — 0.255 for 
Model 1 and g = 0.54 for Model 2). 



where I?o(ripoak) = Max[I?n]. 

Equation (|46p can be used to renormalize a frequency 
landscape numerically, just like an Ising lattice can be 
renormalized by assigning the direction of the block spin 
according to the majority rule. Ideally, 



lim {V^, 



(48) 



but there are numerical problems (apart from the dy- 
namical instability). First, small frequency gradients in 
a presumed cluster are magnified by Eq. (j46p . Second, 
a block oscillator containing a cluster border will not be 
part of the renormalized cluster, which distorts the renor- 
malization of cluster sizes for small and medium sized 
clusters, like those obtained in a simulation. Basically, 
the problem is that the frequencies are continuous vari- 
ables, whereas spins are discrete. 

FigurefTUlshows numerical estimations of for Model 
1 and Model 2. Panels (a) show that fipcak ~ E[fife], at 
least for g < gd, indicating that Eq. ([47)1 is fulfilled. 
This is true even if is not symmetric about its mean. 
However, in Model 2 for g > gd, it is clear that (fipcak) > 



Panels (b) in Fig. [TOlshow the low frequency tails for 
each model at different values of g. Just below and at 
g — 5ci it seems that 



(49) 



with X w 4 (dashed lines), suggesting that the require- 
ment that Var[f2fe] exists at g — gd is indeed fulfilled. 
Regarding the high frequency tails, in Model 1 they seem 
to fall off quicker than a power law for the largest frequen- 
cies for all g. In model 2, there seems to be no tails close 
to and above gd- 

Figure [Til shows numerical renormalizations in phases 
1 and 2 using Eq. The outcomes are qualitatively 

different, indicating that different fixed points are ap- 
proached in the two cases. We have not been able to 
obtain convergence towards the presumed critical fixed 
point density function using a system close to g = gd ■ 
Instead, the outcome is similar to that shown in phase 
1, although the normal distribution is approached more 
slowly (as b increases). The reason may be that, numeri- 
cally, some positive correlations are still remaining (Fig. 

The frequency correlation fimction Tii{r) drops close 
to zero in phase 2, but it seems that it never becomes 
negative. This means that Var[17j] > Var[rife]. Thus the 
fiow under Ri, cannot go towards perfect synchronization 
(r — 1), for which Var[r2fc] — 0, but it can reach states 
where the lattice is synchronized except for isolated oscil- 
lators with opposing frequencies. Such states are indeed 
seen at coupling strengths slightly larger than gc2 (Fig. 
[7]). If there is a non-trivial fixed point corresponding 
to such a state, then the outlier oscillators must have 
a fractal spatial distribution. Otherwise they will "eat" 
the synhronized part of the lattice, and r — ^ as 5 — > oo. 
This effect is seen in Fig. [11] as a (slightly) decreasing 
height of the spike and an elevated baseline of outlier 
oscillators as b increases. 



D. Cluster frequencies 

Assume that the frequency of a cluster C with spatial 
size >> 1 is bounded by the inequality 



\^c - rriool < Ar2i„ax(S'), 



(50) 



where ftk = whenever k ^ C [21]. For t >> 1, choose 
b « S and apply pb- We get S = b^'^S, and fig, — 
ruoo = b^/'^{^c — rrioo) from Eq. ([^ . Consequently, 
An'^^^{b-°S) = 6^/2^f2max(5'), and at a fixed point we 
get 



(51) 



Thus, the cluster frequencies vary less and less as their 
sizes increase. 

Figure [T^ shows comparisons between the theoretical 
prediction in Eq. (|5ip and numerical data. The numer- 



ical difference Q„ 



f^min as a function of S is studied 
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FIG. 12: nmax and nmin are the minimum and maximum fre- 
quencies of clusters of size S found numerically. Let us write 
An = f^max ~ f^min- Lowcr bouuds ou Anniax(5') are shown, 
assuming dAilmax/dS' < 0, VS (picewise linear, continuous 
curves). Dashed lines: predictions by Eq. (I51|l for a critical 
ensemble. Data from 10 realizations of {cuk} for each g. 



rather than the difTerences f^max — Woo or rhao — ^min, 
where rhoo is a numerical estimation of TOoo • The reason 
is that I want to estimate as few quantities as possible. 
Especially for large S, the latter differences are very sen- 
sitive to the choice of rhoo- Note however, that if these 
differences are used, data is obtained that support Eq. 
([STjl to the same extent as the data shown in Fig. [12] 
does. Around gd and in phase 2, the agreement with 
theory is reasonable, given the fluctuations in the data. 
This gives further support to the idea that the entire 
phase 2 is critical. For Model 2, the large values of Ail 
for g = 0.50 and 0.54 for S" = 1 and S = 2 are due to 
the existence of isolated oscillators which are transiently 
suppressed, with Qk ~ 



E. Frequency transient 

From Eq. 1^ we get m{t') = E[d^j/dt'] = b^/'^m{t)- 
(6^/2 „ \)moo- At a fixed point we have m*{h-^/'^-'^t) - 
m*^ = b^/'^{m*{t) — m*^) with solution 



m*{t) 



(52) 



Close to the critical couplings gd ~ 0.23 (Model 1) 
and gd « 0.50 (Model 2), the agreement with Eq. ([52]) 
is excellent in the data shown in Fig. [T31 In the double- 
logarithmic plots, there is a tendency to a gradual in- 
crease of the slope as g increases, suggesting that the 
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FIG. 13: Transient of mean frequency m{t). Lattice size 
2000 X 2000. The dashed lines correspond to the prediction 
by Eq. (|52p for a critical ensemble. Good data quality en- 
ables use of dm/dt, so that no asymptotic values have to be 
estimated (C.f Fig. [T4)) . 



scaling expressed in Eq. (|52p only applies at gd , and not 
in the entire phase 2. This in turn suggests that phase 
2 is not critical, at least that it is not attracted to the 
critical fixed point £* described in this paper. 



F. Mean frequency for finite A*' 



Taking the ensemble mean of Eq. ([36|l gives E[(I'j] 



6^/2 (E[r, 



We may write TOoo = {N) , 



Too 

and then have E[mj] = moc{b^)- If we first let the size 
of the whole lattice go to infinity and then set = N, 
we get 

moo{N) - moo(oo) = N'^^^ {Epj] - moo(oo)} . (53) 

Taking the ensemble mean of Eq. in the limits N — > 
oo and t oo, we get moo(oo) = Ep^] + E[/),^]. Using 



Eq. we may therefore write 



(54) 



At a critical fixed point £* , E*[hj] — E*[/ij,] for all b 
(or N) and therefore moo{N) — moo(oo) oc 

TV- 1/2 for 

all N. At a critical ensemble attracted to £* , we have 
E*[£,] ^ E*[hf,] as & c», so that 



moo{N) - moo(oo) cx TV^^/^, N » 1. 



(55) 



For odd coupling [Eq. we have E[/i^] = E[/ij.] = 

and moo{N) ~ E[a;fe] for all N . This corresponds to zero 
constant of proportionality in Eq. I|55p . 

For sub-critical ensembles it is expected that 
limh^oo E[/i ■] = 0, and for super-critical ensembles that 
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FIG. 14: Mean frequency m(J) at time t = 1000 versus 
L = a/]V. Averages of up to 400000/Af realizations. Esti- 
mated values of m(oo) (dashed horizontal lines) provide good 
agreement with Eq. (|55|l for all g (dashed diagonal lines). 



linih_>oo E[/ij] =00. In both cases, Eq. ([55)1 cannot be 
expected to hold as ^ 00. 

Nevertheless, the data in Fig. [T3] is consistent with 
Eq. (|55p for all shown g. However, crossover to an- 
other scaling for larger L — ^/N cannot be excluded. 
For both models, the asymptotic behavior is reached for 
larger L for higher values of g. Due to poor data quality 
(C.f. Fig. [T31), estimated values of m{oo) have to be re- 
lied upon, chosen to get curves in the double-logarithmic 
plots that are as straight as possible for large L. An- 
other possible source of error is that I had to compute 
the mean frequencies at a rather small time t = 1000 due 
to limited computational resources. However, tests with 
smaller and larger t indicate that this is not crucial. 



G. Correlations of interactions 

Let us turn to the renormalization of the interactions 
hk- It turns out to be useful to decompose hk into the 
coupling functions fik , and to define the asymmetry func- 
tion 



dmix,y) = fikiy,x) + fki{x,y), 



(56) 



where m is the edge connecting oscillators / G j and 
k £ j. Let n be a directed edge across 6j (Fig. [15]) and 
let us write /„ = /a-, where I ^ j and k £ j. We then 
have 



i;;,)], (57) 
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FIG. 15: The interaction hk is a sum of the 2D coupling 
functions fik- In the same way, the interaction hj of block j 
can be decomposed into a sum of 2Db^~^ border terms /„ 
and Db^{l — b~^) interior tems dm ~ fik + fki- 



where d;}^{x,y) = f{i{x,y) + fkiiv^x). 
lowing expressions more compact, let 



Ad 



■^7n —m —m 



To make the fol- 



(58) 



be the mean increase of dm as block j is connected to its 
neighbor block oscillators. For odd coupling, i.e. 



fik{y,x) = ~fki{x,y), yik, 



(59) 



we have dm = dm = 0. An example is the Kuramoto 
model fik{x,y) = sin[27r(x - y)]. 

Information about critical behavior can be gained by 
comparing moments of the original and renormalized in- 
teractions: E[/ij,], E[ft,^], and so on. A comparison be- 
tween E[h^] and E[hj] just leads us back to Eq. (|54|) . Be- 
low, I focus instead on Var[/ij.] and Var[/ij], from which 
information can be gained of two-point correlations of 
the interaction. At this point we make use of assumption 
(|23| . We may then write 



(60) 



Var[^^.] = b-^ J2 ^OY[h,~h^,h,,~h^,], 

k.k'ej 

or, upon decomposition, 

Var[^^.] = b-^g^j:n,n'Cov[l^JJ + 

b'''9'j:n,n.C0v[l^,Adm] + 

^-VE™.™'Cov[Arf„,Ad„,] 
= S1 + S2 + S3. 



(61) 



In the following, three correlation functions will be used: 



r/(0 

rAd(r) 



= hmt^c 



It ^00 Var[/^^j 

C°v[/,..Ad„]|„_„|=,. 
Cov[/^,A<i^)|„_„l = i 
Cov[Ad^,Ad^]|„„„/ 



(62) 



limt- 



Var[Ad^ 



Note that fik has a direction of influence I k, and that 
dm is vertical or horizontal. Correlations between dif- 
ferent types of pairs should therefore be separated, and 
summed up to yield the covariances in Eq. (j6ip . Numer- 
ically, only parallel pairs are considered. 
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Si S2 S3 

"■ft ' 

r 



FIG. 16: Pairs of interactions and distances in a block- 
oscillator j, used in the expressions for Si, S2, and S3 in 
Eqs. dSI}, dSl, dSB, and i^. 



A necessary fixed point condition is Var*[/ij] = 
Var*[/ij.] for any & as t —> 00, or, in particular, at a non- 
trivial fixed point. 



Model 1 



Model 2 



lim lim Var*[/i ] = 
We may write 



const. > 0. 



5*1 oc 6 



-D 



Ni{r)Tf{r)dr, 



where 



iVi(r) ^©(fo^^V^-^) 



(63) 



(64) 



(65) 



is the number of pairs nn' at distance r (Fig. I16|) . It 
follows that 



lim S'l = const. > F ^ (r) cx r 



2-D 



Similarly, 



where 



5*2 oc b 



-D 



N2[r)VjAd{r)dr, 



N2{r) =C'(6^-V^~^) 



(66) 



(67) 



(68) 



is the number of pairs nm at distance r (Fig. I16[) . There- 
fore it is expected that 



l-D 



lim 5*2 = const. > Fi^^(r) cx r 

Turning to 6*3, we may write 

Cov[Ad„,Ad„,] = V(/o)rAd(r), 



(69) 



(70) 



where p is the smallest distance from any of the two edges 
m or ml to (5j, and r is the distance between m and m! 
(Fig. [TH). The function 



-0(p) = Var[Ad„ 



(71) 



measures how the lattice that surrounds block j, in the 
mean, changes at distance p from (5j, so that we have 



V'(l) > 



hm,. 



>^(p) -0. 



(72) 



3 



-10 



g = 0.20 - 



-6.5 



ln(p) 




ln(p) 



FIG. 17; ^(p) = Var[Ad„] as a function of the distance p 
from 8j to the edge m (Fig. [16}. For both models, it seems 
that ■!/'(p) (X p~" at gd. For Model 1, a ^ 1/4, and for Model 
2, a « 1/2. Lattice size: 140 x 140. Block size: 70 x 70. 
Transient time: t = 5000. For each of 10 realizations of {uik}, 
20 initial conditions 0(0) were used. 



We may then write 

/b/2 pb-2p 
ij{p) J N3{p,r)TAd{r)drdp. (73) 



Here, 



^3(p,r)-0(p^-V^-i) 



(74) 



is the number of pairs mm' for given p. In a critical fixed 
point ensemble, it is expected that 



V'*(p)«p-" 

and therefore we get 

lim 53 = const. > <^ F 



b — *oo 



Ad 



(r) cx r 



a-D 



(75) 
(76) 



provided a 1. 

I have not been able to deduce the value of a from 
first principles. Figure [17] shows numerical estimations 
of tp{r) for Model 1 and Model 2. It seems that a — 1/4 
for Model 1, and a = 1/2 for Model 2, and thus that 
it is a non-universal, model dependent critical exponent. 
The scaling form ([75]) seems to apply only at 17 = gci, 
suggesting that phase 2 is not critical. 

It can be argued that since d^ and Ad^ are linear com- 
binations of /^^ and f^^, Tf, F^, F/Ad and Fa^, should 
have the same functional form for r >> 1. Then, 



F^cxFScxr-'^ 



with 



l3 = D-2, odd fik 

(3 = D - I, other fik, a>l 

(3 = D — a, other fik, a < 1 



{case 1) 
(case 2) 
{case 3). 



(77) 



(78) 



The terms Si , S2 and ^3 are responsible for criticality in 
the three cases, respectively. In case 1, 



SI ~ const. > 



^2,3 - 0. 



(79) 
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g = 0.20 
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0.54 




0-50A • 




g = 0.45 , m 







ln(r) 



Model 1 





-0.2 L 



^ Vf-. g = 0.30 

5: g = 0.23 

r^: g = 0.23^'" 







100 



FIG. 19: Illustration of the fact that numerical estimations 
of Tf{r) are less well-behaved than those of rd(r). For a 
given finite lattice size, F / (r) becomes more ill-behaved when 
g increases (see text). In Model 2, a time average of each fik 
was calculated during At = 100. (C.f. Fig. [T8)l . 



FIG. 18: Pair-correlations of / (left) and d (right). According 
to Eq. (|78|) and the estimations of a in Fig. 1171 it is expected 
that (3 = -7/4 for Model 1, and (3 = -3/2 for Model 2 in 
critical ensembles (dashed lines). Ten 4>{0) were used for each 
g to estimate the initial condition mean. In Model 2, a time 
average of each fik was calculated during At — 100, due to its 
pulse-like nature. The choice of At did not affect the results. 

[In fact, 52,3 = for all ensembles since djn{x) = 0.] In 
case 2, 

limb„,oo S2 = const. > , , 

linib^oo 3 = 0, ^ ^ 

and in case 3, 

limb^oo 5"! = const. > , . 

limb^oo SI 2 = 0. ^ ' 

In cases 1 and 2, critical behavior is ruled out below 
D — S and D = 2, respectively, since correlations must 
decay with r. In case 1, the result -De > 2 by Daido [2^ 
is regained. 

The numerical results in Fig. [T7] suggest that a < 1 
for both Model 1 and Model 2, and thus that the term 
53 is responsible for criticality in both models. However, 
since the estimated values of a differ, it is possible that 
there are other non-odd models whith a > 1, in which 
case ^2 becomes the crucial term. 

Figure [TH] shows numerical estimations of F^ (r) and 
Td{r). Looking at Td{r), the data is consistent with the 
combined theoretical and numerical predictions [Eq. (|78p 
and Fig. [TT] in both Model 1 and Model 2. Looking at 
F/ (r) , the data are consistent with theory only in Model 
27 

The reason for this discrepancy between numerical 
data and theory in Model 1 is likely to be found in the 
fact that F/(r) is less well-behaved than Fd(r) for finite 
lattice sizes. The phase fields 4'{x, y) become more and 
more well-ordered as g increases, containing just a few 



foci or spirals as phase 3 is approached (at the present 
lattice size) ^24] . Therefore the phase waves tend to move 
in opposite directions at opposite ends of the lattice, giv- 
ing rise to negative correlations of / at large distances. 
This dependence on the wave direction of the correlations 
is eliminated by the definition of d [Eq. (|56p]. 

This problem is illustrated for Model 1 in Fig. [121 
Close to criticality, for g = 0.23, Fd(r) converges nicely 
towards zero as r increases, whereas Yf{r) drops signif- 
icantly below zero, and then fluctuate, at least up to 
r = 150. (This is the maximum r considered, since the 
lattice size is 300 x 300.) This effect is more prominent for 
larger g as seen in the estimation of Fy(r) for g = 0.30. 
The zero-crossings of F^ (r) is the reason why the curves 
drop sharply in the double-logarithmic plots in Fig. [T8l 



H. The correlation length 

Let us analyze F/ close to a critical fixed point in a 
subcritical ensemble, and make the standard ansatz 

F^(r) = cr-'3e-''/«('^s)^ (82) 

where 

Ag = (g-g*)/5*. (83) 

As discussed below, subcriticality is expected only for 
g < g* . It is therefore assumed that Ag < 0. We may 
write Var[/j^] = F{g)- Assuming that dF/dg 7^ at 
g = g* , we have 

AgcxVar*[/J-Var[/J (84) 

for small enough Ag. Consider the case of odd coupling. 
From the expression for Si we get 

Ag'oib-^J r^-^[T}{r)-Tf_ir)]dr. (85) 
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Taylor expanding the exponential part of F/ gives Aji' oc 
b£,~^. Using ^' = ^/5, specifying £^ — it is seen that 
the correlation length of the initial condition mean of the 
coupling fik diverges according to 

C/cxAg-i, (86) 

for small enough \Ag\ if Ag < 0. A similar calculation 
gives the same result in the cases of non-odd coupling, 
using the expressions for S2 and 5*3. 

Indeed, simulations suggest that ^/ diverges at 5 = gd, 
but unfortunately I have not been able to obtain good 
enough data to test relation The fluctuations in 

the estimated ^/ are too large close to gd- (I used up 
to three realizations of {ujk} for each g, and for each 
{wfc}, ten (j)(0) were used to estimate the initial condition 
mean). It was not possible to use data from estimations 
of Td either, since it drops close to zero for too small r 
to be able to resolve its functional form. 



I. Direction of the renormalization flow 

In Fig. [17] the exponent a is estimated in the rela- 
tion ip{p) cx [Eqs. (|7T|) and (|75| ]. that is expected 
to hold in a critical ensemble. Let us call these estima- 
tions ai and a2 for models 1 and 2, respectively. In 
phases 2 and 3, ip{p) clearly falls off slower than this. 
Figure [18] shows that in phases 2 and 3, F/ and F/ 

falls off as r-(^-"i) (Model 1), r-'-"'"^^ (Model 2), or 
possibly slower. Taken together, these observations sug- 
gest that condition ([76]) is violated in phases 2 and 3, 
and that limb^oo S3 — 00. This in turn means that 
linih^oo Var[ft,j] limh_>oo Var[ft,'j,] = 00, and that the 
renormalization flow goes in the direction of increasing g 
for g > gd (Fig. [20]) . That the flow goes towards g — 
for g < gc becomes clear from a similar argument. I have 
mentioned the possibility that the entire phase 2 is crit- 
ical, and that it is attracted to the critical fixed point 
£* . Some numerical results favor such an interpretation 
(see Figs. [3 [H [H [H [H and also Ref. H^). How- 
ever, based on the combined numerical and theoretical 
argument given above, I hypothesize that this is not so. 

Referring to the discussion in section IV C[ it seems 
that the renormalization flow in phase 2 cannot approach 
states with r = 1. Therefore it is probable that phase 2 
is invariant under i?^. There may be a second, attractive 
fixed point with Var[f2i:] = 00 somewhere along the line 
separating phases 2 and 3, possibly at infinity where g — > 
00 or CT^ ^ 00. 

Correlation functions seem to decay as a power law 
or slower for all g > gd- In fact, Eqs. ([M]l . ([57]) and 
(170)) predict that finite correlations lengths are excluded 
for g > gd, since whenever Tf has an exponential factor, 
limh^oo Si-3 = 0. This corresponds to a renormalization 
fiow towards g = (Var[^^] — 0), which can be expected 
only for g < gd- Therefore, phases 2 and 3 must be 
considered supercritical. 



< + ^ \ > -9 

9c^ 9c2 
Phase 1 Phase 2 Phase 3 

(subcritical) (supercritical?) (supercritical) 



FIG. 20: Hypothetical flow under Ri, in a one- dimensional 
projection of {£}- The critical coupling strength gci is sup- 
posed to belong to the stable manifold Sc of the critical flxed 
point £* (c.f. Fig. [2}. It is speculated that the flow does not 
pass gc2, i.e. that phase 2 is invariant. See text for explana- 
tion. 



VI. DISCUSSION 

In this paper, I present a real-space renormalization 
transformation for oscillator lattices with quenched disor- 
der. The transformation acts on ensembles of lattices and 
predicts the behavior of ensemble averaged quantities. It 
is assumed that the variance of the natural and attained 
frequencies exists, but it should be possible to generalize 
the theory. A bold hypothesis is that if a system of form 
Q is critical for some parameter values, then the critical 
behavior is given by the critical fixed point £* desribed 
in this paper. At its present stage, the theory cannot be 
used to decide whether a given system posesses a criti- 
cal phase transition. However, lower bounds on critical 
dimensions for different classes of systems are given. 

In this respect, the crucial difference between odd and 
non-odd coupling stands out clearly in the analysis. For 
non-odd coupling, macroscopic synchronization cannot 
be ruled out for any dimension D > 1, whereas for odd 
coupling it is necessary that D > 3. Perfectly odd cou- 
pling must be regarded as a non-generic special case, ex- 
cept for particular problems that can be mapped onto 
Kuramoto-like models, such as Josephson junction arrays 

The merits of the approach are that it is simple, that 
it applies to a broad class of systems, that several predic- 
tions about critical behavior can be extracted, and that 
it is potentially exact. Most of the predictions have been 
tested numerically with two structurally different two- 
dimensional models. The agreement with theory ranges 
from acceptable to very good. The drawback of the ap- 
proach is that the theory must be considered heuristic at 
its present stage. Its full potential and its mathematical 
foundation should be clarified. 

My experience is that it is computationally demand- 
ing to get good numerical data to compare with theory. 
Large oscillator lattices [0(10^) oscillators] and long sim- 
ulation times [O(IO^) periods of oscillation] are typically 
needed to see critical behavior. Further, to get good en- 
semble averages, it seems that 0(10) realizations of the 
initial condition are needed for each of 0(10) to O(IOO) 
realizations of natural periods. In other words, O(IOO) 
to 0(1000) realizations are needed for lattice sizes and 
integration times of the above order og magnitude. This 
is probably the reason why almost no clear-cut numerical 
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results regarding the existence or non-existence of phase 
transitions in oscillator lattices have been presented in 
the past (section IIIB|) . The data presented in this paper 
should be seen as an initial overview of the behavior of 
some relevant quantities. A more detailed study of each 
quantity is needed. In particular, the number of realiza- 
tions of natural periods has to be increased. 

To put the theory to further test, it goes without say- 
ing that simulations of oscillator lattices with dimensions 
other than D = 2 are called for. Perhaps the quanti- 
ties used in this paper can be used to find an answer to 
the long standing question whether there is a transition 
to macroscopic synchronization in the three-dimensional 
Kuramoto model. 

It is worth noting that the relevance of the second crit- 
ical coupling gc2 is established in this study. It was first 
described in Ref. [1^ , but there a density function of 
natural frequencies with finite support was used. Here, I 
find that it is present even if has tails. It is therefore 
a more generic transition than that to i? = 1 in the glob- 
ally coupled Kuramoto or Winfree models [131, appear- 
ing when has no tails. The nature of the transition at 
gc2 is a subject for future work, and the question whether 
there is an additional non-trivial fixed point associated 
with this transition is left unanswered. 

Theory and numerics taken together indicate that the 
renormalization flow goes towards increasing g for g > gd 
(Fig. [20]) . I judge that both phase 2 and phase 3 are su- 
percritical, and in section IVIl I give a technical argument 
why this is so. Here, a qualitative argument is presented 
why an oscillator lattice cannot be subcritical above gd , 
that is, why correlation functions cannot have exponen- 
tial tails, corresponding to finite correlation lenghts. 

Correlation lengths relate to the typical distance a per- 
turbation or fluctuation spreads. Let us compare with 
the Ising model, which is subcritical both above and be- 



low the critical temperature Tc. In the ordered phase 
below Tc, most spins are aligned. Let us introduce a per- 
turbation in the form of a spin with opposite direction. 
Such a spin increases the probability that a neighbor spin 
will also flip. The perturbation tends to spread. How- 
ever, the lower the temperature, the smaller the proba- 
bility that the neighbor will flip, according to the Boltz- 
mann distribution. Thus, a typical perturbation spreads 
shorter distances, and the correlation length drops. 

The situation is quite different in the ordered phase of 
an oscillator lattice, where I am thinking mainly of states 
with partial frequency synchronization (0 < r < 1). A 
perturbation in such a lattice corresponds to an oscilla- 
tor k that runs at a different frequency. This perturba- 
tion spreads to the rest of the lattice via the coupling 
functions, which will not vary with the entrained fre- 
quency. Assume for simplicity that the coupling has the 
form gipiki<t>i — 4>k)- The peak magnitude of this pertur- 
bation can only increase with g, since the argument takes 
on all values in the range [0, 1) because fi; and fi^ are 
assumed to be different. Thus, if the correlation lengths 
are infinite at a critical coupling gd, they should stay 
infinite even ii g > gd- 

In conclusion, I hope that this study will inspire fur- 
ther theoretical and numerical work on macroscopic syn- 
chronization in oscillator lattices. Unfortunately, the un- 
derstanding of these systems has fallen way behind the 
understanding of globally coupled oscillator networks. 
A better understanding of oscillator lattices should also 
promote the understanding of transitions to macroscopic 
synchronization in complex networks, since the topology 
of these often can be seen as lying in between the topolo- 
gies of the lattice and the fully connected network. 

This study was supported in part by the Royal Swedish 
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